function [dxi] = rcnl_dxi( ...
    dparams,            ...
    s_jt,               ...
    xlinear,            ...
    xnonlin,            ...
    cdindex,            ...
    dfull,              ...
    IVD,                ...
    W,                  ...
    IVD_alt,            ...
    W_alt,              ...
    fe,                 ...
    expmval_mat,        ...
    rho_transform,      ...
    tolerance           ...
)

Pi  = dparams(1:end-1);
rho = dparams(end);
if rho_transform
    rho = rho_transform * (1/(1+exp(-rho)));
end

[delta, ~] = rcnl_meanval(Pi, rho, s_jt, xnonlin, cdindex, dfull, expmval_mat, 0, tolerance);
ddelta = rcnl_ddelta( ...
    Pi,            ...
    rho,           ...
    delta,         ...
    xnonlin,       ...
    cdindex,       ...
    dfull,         ...
    rho_transform  ...
);

[~, dxi] = rcnl_jacobian(ddelta, xlinear, IVD, W_alt, fe, IVD_alt);
